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Abstract 

We have determined the ground-state energies of para-H2 clusters at zero temperature using 
the diffusion Monte Carlo method. The liquid or solid character of each cluster is investigated by 
restricting the phase through the use of proper importance sampling. Our results show inhomoge- 
neous crystallization of clusters, with alternating behavior between liquid and solid phases up to 
N = 55. From there on, all clusters are solid. The ground-state energies in the range N = 13-75 
are established and the stable phase of each cluster is determined. In spite of the small differences 
observed between the energy of liquid and solid clusters, the corresponding density profiles are sig- 
nificantly different, feature that can help to solve ambiguities in the determination of the specific 
phase of H2 clusters. 
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I. INTRODUCTION 



Molecular para-hydrogen has been suggested from long time ago as the possible second 
natural superfluid after liquid 4 He.- However, the H2-H2 interaction is so deeply attractive 
that it crystallizes before arriving to the superfluid transition temperature. The advantage 
of having half the mass of 4 He is therefore not enough to compensate the hydrogen strong 
attraction and para-H2 becomes solid at temperature T = 13.96 K. Many attempts to 
supercool liquid hydrogen down to the expected lambda transition (T\ = 1 — 2 K) have been 
unfruitful, at least for the bulk phase.- Partial success has only been achieved in confined 
geometries, mainly in small pure liquid drops 3 or bigger drops in a 4 He environment. 4 

Para-H2 clusters have been the object of many studies in the last years**""— due to the 
primary interest of determining i) its liquid or solid character and ii) the dependence of 
superfluidity on the size and temperature of the cluster. The first path integral Monte Carlo 
(PIMC) simulation of H 2 clusters, carried out by Sindzingre et al.,- showed that clusters 
comprising up to a number of molecules iV ~ 18 are superfluid at temperatures below T = 2 
K. The maximum number of molecules that shows superfluid behavior has been enlarged 
up to iV ~ 26 in a recent PIMC simulation where lower temperatures T = 0.5 K have been 
analyzed.- The results obtained in Ref.- show evidence of superfluidity mostly localized in 
the surface of the cluster, pointing to an inhomogeneous structure with an inner solid core 
surrounded by a liquid skin that, at the lowest temperatures, is superfluid. This localized 
superfluidity has been questioned in a recent PIMC calculation where it has been shown that 
superfluidity is a global property of the cluster in spite of its significant spatial structure.— 
In the limit of zero temperature, the structure and energy of small H 2 clusters have been 
accurately studied using both diffusion Monte Carlo (DMC)- 1 ^ and path integral ground 
state (PIGS). 7 Both at finite and zero temperature the simulations show the presence of 
magic-cluster sizes in which the chemical potential shows a kink. These more stable N 
configurations have also been observed experimentally in cryogenic free jet expansions using 
Raman spectroscopy. 3 

In the present work, we address the question of the solidification of para-H 2 clusters in the 
limit of zero temperature and as a function of the number of molecules. Our aim is the screen- 
ing of the more stable ground-state structure by performing simulations where the phase 
(liquid or solid) is kept fixed. To this end, we use the diffusion Monte Carlo (DMC) method 



that solves stochastically the iV-body Schrodinger equation exactly for bosons, within some 
statistical errors. Differently from previous studies using DMC, we have focused our atten- 
tion on the discrimination between liquid and solid phases as a function of the number of 
molecules. This goal is achieved by carrying out parallel simulations for liquid and solid 
configurations at each N, with special effort on the search for optimal lattices on which to 
build trial wave functions for the solid clusters. With the present study, we show which is 
the energetically preferred phase at T = for each N, in the range 13 < N < 75, and how 
the energy difference between the two type of clusters changes with increasing N. 

The rest of the paper is organized as follows. In the next Section, we present our theo- 
retical approach that relies on the DMC method. In Sec. Ill, we report our results on the 
energetic and structure properties of the para-H2 clusters in the N range analyzed. Finally, 
we end with the summary and main conclusions in Sec. IV. 

II. QUANTUM MONTE CARLO APPROACH 

Our study of para-H 2 clusters relies on a purely microscopic approach whose inputs are 
only the interparticle interaction and the mass. Our goal is the study of these finite systems 
at zero temperature to deal with their ground state. To this end, we use the DMC method 
which is able to generate exact (within statistical uncertainty) information through guided 
random walks. The starting point in the DMC method is the iV-body Schrodinger equation, 
written in imaginary time 



with R = (ri, . . . , rjv), a 3iV-dimensional vector (walker), and t the imaginary time measured 
in units of h. The time-dependent wave function of the system \I/(R, t) can be expanded in 
terms of a complete set of eigenf unctions 0i(R) of the Hamiltonian, 



where E^ is the eigenvalue associated to </>j(R). Consequently, the asymptotic solution of [U 
for any value E close to the energy of the ground state and for long times (t — > oo), gives 
0o (R), provided that there is a nonzero overlap between \1/(R, t = 0) and the ground-state 
wave function 0o(R)- 




(1) 




(2) 



n 



A direct Monte Carlo implementation of [T] is hardly able to work efficiently, especially 
when the interatomic potential contains a hard core. This is solved by using importance 
sampling. The importance sampling technique is a general concept in Monte Carlo and is 
one of the best methods to reduce the variance of any MC calculation. The importance 
sampling method, applied to [TJ consists in rewriting the Schrodinger equation in terms of 
the wave function 

/(R,*) = V(R)*(R,*), (3) 

ip(R) being a time-independent variational wave function that describes approximately the 
ground state of the system. Considering a Hamiltonian of the form 



[U turns out to be 
df(R,t) 



-D V R / (R, t) + D V R (F(R) /(R, t) ) + (E L (R) - E) /(R, t) , (5) 



dt 

where D = h 2 /(2m), E L (R) = ijj(R)~ 1 Hil)(R) is the local energy, and 

F(R) = 2^(R)- 1 V R ^(R) (6) 

is called drift or quantum force. F(R) acts as an external force which guides the diffusion 
process, involved by the first term in[5l to regions where t/*(R) is large. 

The r.h.s. of [5] may be written as the action of three operators Ai acting on the wave 
function /(R, t), 

df(R,t) 



dt 



(A 1 + A 2 + A 3 )f(R,t)=Af(R,t) (7) 



The three terms A4 may be interpreted by similarity with classical differential equations. 
The first one, A\, corresponds to a free diffusion with a a diffusion coefficient D; A 2 acts as 
a driving force due to an external potential, and finally A 3 looks like a birth/death term. 
In Monte Carlo, the Schrodinger equation [7] is best suited when it is written in an integral 
form by introducing the Green function G(R',R,t), which gives the transition probability 
from an initial state R to a final one R' during a time t, 

f(R', t + At) = J G(R, R, At) f(R, t) dR . (8) 



More explicitly, the Green function is given in terms of the operator A = A% + A2 + A3 by 

G(R,R, At) = (R'| exp(-AAt) |R). (9) 

DMC algorithms rely on reasonable approximations of G(R', R, At) for small values of 
the time-step At. We work with a second-order expansion of the exponential [9] to reduce the 
time-step dependence.— Once a short-time approximation is assumed, M is iterated repeat- 
edly until reaching the asymptotic regime /(R, t — > 00), a limit in which one is effectively 
sampling the ground state. 

Para-H2 is a boson particle with total spin and with rotational ground-state state 
J = 0. It is well described by a merely radial interaction due to its high-degree sphericity. 
Among the different interatomic potentials proposed for describing the H 2 -H 2 intermolecular 
interaction, we have chosen the Silvera- Goldman potential^ due to its proved accuracy and 
its dominant use in microscopic calculations as the present one. The well depth of the 
molecular hydrogen interaction is ~ —37 K, a factor of four larger than in helium, making 
the ground-state phase of bulk at zero temperature be an hep crystal in spite of H 2 mass 
being half the 4 He one. 

The trial wave function used for importance sampling is written as the product of one- 
and two-body correlation factors, 

N N 

r/)(R) = J] e u{r '^ \{e v{r ^ , (10) 

l=i<j i=l 

with 

v (r) = -ar 2 . (12) 

The two-body term u{r) accounts for the correlations induced by the potential V(r) and 
also for the finite size of the system that implies the wave function approaches zero when 
the distance is of the order of the cluster size. The one-body term v(r) is only used for solid 
clusters and localizes particles around the preferred sites (capital indexes in fTUi) . This model 
wave function for the solid, with the one-body Gaussian terms, is the well-know Nosanow- 
Jastrow (NJ) wave function that has been widely used in the study of bulk quantum solids. 
The NJ wave function is not symmetric under the exchange of particles but the influence 
in the energy of symmetrization is known to be of the order of miliKelvin^ and therefore 



indistinguishable within our statistical errors. Recently, it has been proposed a symmetrized 
NJ wave function to study superfluidity in bulk solid 4 He and proved that the change in 
energy due to Bose symmetry on top of the NJ model is imperceptible, within the numerical 
resolution of quantum Monte Carlo.— It is worth noticing that the effect of symmetrization 
in the structure of para-H 2 clusters at temperature T = 1.5 K has been recently studied by 
Warnecke et a/.;— the results obtained show a very small influence of Bose statistics on the 
density profiles, the energy differences being not reported. 

The parameters entering in Eq. (11) and Eq. (12) have been optimized using the varia- 
tional Monte Carlo (VMC) method. For liquid configurations the optimal parameters are: 
b = 3.58 A, = 2.79 A -1 , and a = 0; for solid ones: b = 3.32 A, fi = 0, and a = 0.521 A" 2 . 
These are the VMC optimal parameters for N = 40 but the dependence of these parameters 
with N is tiny: for N = 13, = 2.30 A -1 , a = 0.360 A~ 2 , and b is the same. In DMC we 
have neglected this slight N dependence because the results are insensitive to it. On the 
other hand, two technical issues related to the implementation of the DMC method, i.e., the 
time-step dependence and the number of walkers, have been accurately analyzed to reduce 
any systematic bias to the level of the statistical noise. 

The simulation of solid clusters with the NJ wave function we are using requires of a set 
of lattice points. In the case of bulk solids, the number of possible lattices is few, well known 
and easy to characterize geometrically. When dealing with finite systems, the issue of the 
best geometrical arrangement is somehow universal for the smallest structures, where Mackay 
polyhedra are the preferred structures, but it becomes more complex when the number of 
particle increases. Moreover, the determination of the optimal solid patterns are normally 
obtained using classical physics and the influence on those of quantum derealization is less 
known. Our strategy for this optimization has been the use of both simulated annealing (SA) 
and ab initio random search (method for finding structures where the total classical force 
felt by any particle is approximately zero). The structures of minimum energy predicted for 
both approaches essentially coincide for N < 30, but for larger iV values the SA method 
has proven to be better. In the SA approach, we have used both an exponential annealing 
schedule and another with constant thermodynamic speed.— The search for any N has been 
carried out by starting from 20-40 random configurations and selecting the best ones for a 
posterior quantum simulation. In some cases, we have performed additional optimizations 
using long classical Monte Carlo simulations. Finally, we used the best classical solutions 



as initial setup for the quantum simulations and reported in all cases the lowest energy 
achieved. It is worth mentioning that we introduced a scaling factor / s in the quantum 
simulation, 

r° = (r< - r CM )/ s , (13) 

with r CM the center of mass of the cluster and Tj (i = 1, . . . , N) the SA points, to allow 
for possible spatial expansions or contractions of the classic lattice points but the optimal 
energy was always the corresponding to the classic solutions (/ s = 1). 

III. RESULTS 

We have calculated the energy and structure properties of H 2 clusters in the range 13 < 
iV < 75 using the DMC method discussed in the preceding Section. At each N, we have 
used both the liquid and solid trial wave functions introduced in Sec. II. In[p, we show the 
results for the energy per particle obtained for both phases and for clusters up to N = 55 
molecules. The difference between both configurations is small in all the N range studied 
pointing to a highly correlated liquid, even for the lightest clusters. The energy per particle 
in absolute value increases with N for both phases but solid clusters have energies with a 
kind of zigzag dependence whereas the liquid ones follow a smoother law. The irregular 
behavior of the energy of solid clusters is consequence of the appearance of magic numbers, 
with geometrically more compact structures that make them more stable. 

The difference between the energies of liquid and solid clusters is shown explicitly in [21 
Starting from N = 13, the energy of liquid clusters is clearly preferred but the difference 
(El — Eg)/N decreases when N increases. According to our results, the first cluster that 
is solid in its ground state is the one with N = 32. From N = 33 to 40, there are some 
solid clusters but for N = 41 and the following ones there is again a liquid regime. Arriving 
to iV = 55, one enters definitively in a preferred solid phase that persists up to N = 75 
which is the largest cluster here analyzed. The more intriguing aspect of our results is the 
existence of a liquid stability island for N = 41-50 between the initial solid regime and 
the second one, that we think is the stable one for N > 55. It can be argued that the 
nonuniform crystallization of H2 clusters that we have observed can be due to our model of 
solid clusters. We can not discard completely this argument but we have performed a rather 
exhaustive search of lattice points, as commented in the preceding Section, and the results 
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Figure 1. (Color online) Energy per particle of H2 clusters as a function of N. Circles and 
diamonds stand for liquid and solid phases, respectively. The lines are guides to the eye. Error 
bars are smaller than the size of the symbols 

remain unchanged. 

An important part of our simulation has been the search of optimal structures for solid 
clusters. The most useful tool has been the simulated annealing algorithm that, in spite 
of being completely classical, has proven to be able to generate the best configurations. A 
simple check for that has been the optimization of a scaling parameter / s that we introduced 
to make possible contractions or expansions of the classical solution. Systematically, the 
optimal solution has been / s = 1. In|3l the best lattice points for solid clusters with N = 18, 
19, and 20 are shown. We have joined with lines the planar structures to get a better 
visualization of the clusters. The case N = 19 corresponds to a well-know magic number 
and it is in fact a structure that we identify in a lot of clusters. As one can see in [31 in the 
clusters N = 18 and N = 20 the one-defect and one-excess structures, with respect to the 
geometrically perfect N = 19 lattice, are clearly observed. 

The effort of achieving optimal structures increases significantly with N, so a lot of SA 
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Figure 2. (Color online) Difference between the energies per particle of liquid and solid clusters 
((E L - E s )/N) as a function of N. 

simulations with different initial configurations have been carried out. Our best configura- 
tions coincide with reported ones for Lennard- Jones interactions up to iV ~ 30.— For larger 
iV values, our results differ significantly from the published ones^ and we get systematically 
better energies with our structures. In HI we report our optimal lattice points for clusters 
with N = 31 and N = 34, i.e., beyond the regime where our predictions are compatible with 
the ones from Ref.—. As commented before, it is quite remarkable that the inner structure 
of these two clusters is the same and coincides exactly with the magic structure of N = 19. 
Around this well-defined structure it starts to appear a second pentagonal shell which is 
only complete in the center and concentric with the central pentagon of the N = 19 cluster. 

In Table 1, we report the ground-state energies for each N and the corresponding phase. 
Our results for the liquid clusters agree with the previous results of Ref.— obtained also 
using the DMC method and the same interaction potential. The present results for the 
solid clusters are new and never calculated before. The energies reported in Table 1 lead to 
persistence of liquid character up to relatively large N values, in particular to larger values 
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Figure 3. (Color online) Optimal distribution of lattice points for solid clusters with N = 18, 19, 
and 20. 
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Figure 4. (Color online) Optimal distribution of lattice points for solid clusters with N = 31 and 
34. 

than the ones reported in previous PIMC estimations at finite temperature.- 1 ^ Nevertheless, 
the difference in energy between the two phases remains small even for the largest clusters 
studied. We made several attempts at simulating clusters with an inner solid core and a 
liquid surface but we were not able to get any improvement of the energy with respect to 
the optimal values reported in the table. 



Table I. Energy per particle of the ground state of para-H 2 clusters as a function of N. The phase 
of the cluster is labeled L and S for liquid and solid, respectively. Figures in parenthesis are the 
statistical errors. 
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Figure 5. (Color online) Chemical potential of H2 clusters as a function of N. 

The possibility of magic clusters, with higher probability of being experimentally ob- 
served due to its enhanced stability, has deserved the attention of all the theoretical and 
experimental works in this subject. If such structures exist, one can see its signature in 
energy differences, mainly in the chemical potential defined as 

fi(N) = E(N) — E(N — 1) . (14) 

In 0, we show the results for the chemical potential [2] obtained from the DMC ground- 
state energies reported in Table 1. In the liquid regime (small clusters), /i(iV) is quite a 
smooth function with small local minima for N = 23, 25, and 30. For N > 30, a clear 
zigzag structure is observed with some minima that correspond to liquid clusters and other 
that are solid. Therefore, the chemical potential we obtain is no more a smooth function as 
derived in previous DMC estimations where only liquid clusters where considered.- 1 ^ Our 
results show a more complex structure with some alternating phases and with a lot of local 
minima, mainly when entering in the solid phase. 

The issue of the existence of magic clusters can be better analyzed by calculating the 
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Figure 6. (Color online) Second energy difference A2 of H2 clusters as a function of N. 

second energy difference, 

A 2 (iV) = E(N + 1)-2E(N) + E(N-1) . (15) 

DMC results for A 2 (iV) are shown in Fig. |6j The observed pattern emphasizes the behavior 
of n(N), with a quite smooth behavior up to iV ~ 30 and sharp peaks beyond this value. 
By looking at the position of the maxima, we observe magic structures for N = 13, 19, 23, 
25, 28, 30, 34, 37, 40, 43, 47, 49, 51, and 53. Some of these values have been recently found 
in PIMC simulations at finite temperatures.— 1 ^ 

DMC simulations serve also to calculate structure properties of the clusters, as for in- 
stance, the density profiles. In [71 we show the density profiles of clusters with iV = 19 and 
N = 29 using both the liquid and solid trial wave functions. As one can see, the size of the 
cluster is very similar for liquid and solid clusters with the same N but the inner structure 
is rather different. This is particularly clear for N = 19 where the central density of the 
liquid is large whereas for the solid is zero. This central zero density for the solid cluster 
can be well understood by inspection of the lattice points of the N = 19 cluster reported in 
[3j When iV increases, the differences between liquid and solid clusters diminish but are still 
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Figure 7. (Color online) Density profiles of H2 clusters with N = 19 and 29. Labels L and S stand 
for liquid and solid clusters, respectively. 

clear for the case N = 29 reported in [7J It is worth noticing that the density profiles of the 
liquid clusters present a two-shell effect- due to the strong intermolecular H 2 -H 2 attraction 
that is never observed in liquid 4 He clusters. 

We report the evolution of the density profiles of liquid and solid clusters with N in [H] 
and EJ respectively. Increasing N in the liquid clusters produces a continuous change in the 
density profile: the central density decreases to zero and a two-shell structure, that moves 
progressively to larger r, clearly appears. In the case of solid clusters the evolution 
with iV is not so smooth, mainly for the lowest N values. Different from the liquid clusters, 
solid ones show always empty density in the center of the cluster and additional structure 
between the two shells that tend to dissappear when N becomes larger. Comparing density 
profiles for the same N and different phase, one concludes that functions p(r) are in all 
cases different enough to be unambiguously discerned in spite of the fact that the difference 
between their binding energies is less that 1 K per particle. 
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Figure 8. Density profiles of liquid H 2 clusters with N = 19, 21, 23, 25, 27, and 29. The line for 
each N can be identified looking at the radius of the cluster that increases monotonically with N. 



IV. CONCLUSIONS 



H 2 clusters in the range N = 13-75 have been microscopically characterized at zero 
temperature using the DMC method. Controlling the phase of the cluster by using different 
models for trial wave functions used for importance sampling we have been able to determine 
the ground-state stable phase for each N in the range analyzed. Our results point to a 
nonuniform crystallization of the H2 clusters, with some alternating behavior between the 
two phases depending on the particular N value. For clusters with N > 55 the stable phase is 
the solid one. The structure of the clusters, as shown in their density profiles, is significantly 
different for liquid and solid clusters with the same N, even when the difference in energy 
between both is really tiny. Therefore, the shape of the density profiles can help to identify 
the nature of these clusters both in experiment and in finite-temperature simulations. 
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Figure 9. Density profiles of solid H2 clusters with N = 19, 21, 23, 25, 27, and 29. The line for 
each N is of the same type than in Fig. 7. 
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